Strong interactions between highly dynamic lamina-associated domains and the nuclear envelope stabilize the 3D architecture of Drosophila interphase chromatin

Background Interactions among topologically associating domains (TADs), and between the nuclear envelope (NE) and lamina-associated domains (LADs) are expected to shape various aspects of three-dimensional (3D) chromatin structure and dynamics; however, relevant genome-wide experiments that may provide statistically significant conclusions remain difficult. Results We have developed a coarse-grained dynamical model of D. melanogaster nuclei at TAD resolution that explicitly accounts for four distinct epigenetic classes of TADs and LAD–NE interactions. The model is parameterized to reproduce the experimental Hi-C map of the wild type (WT) nuclei; it describes time evolution of the chromatin over the G1 phase of the interphase. The simulations include an ensemble of nuclei, corresponding to the experimentally observed set of several possible mutual arrangements of chromosomal arms. The model is validated against multiple structural features of chromatin from several different experiments not used in model development. Predicted positioning of all LADs at the NE is highly dynamic—the same LAD can attach, detach and move far away from the NE multiple times during interphase. The probabilities of LADs to be in contact with the NE vary by an order of magnitude, despite all having the same affinity to the NE in the model. These probabilities are mostly determined by a highly variable local linear density of LADs along the genome, which also has the same strong effect on the predicted positioning of individual TADs -- higher probability of a TAD to be near NE is largely determined by a higher linear density of LADs surrounding this TAD. The distribution of LADs along the chromosome chains plays a notable role in maintaining a non-random average global structure of chromatin. Relatively high affinity of LADs to the NE in the WT nuclei substantially reduces sensitivity of the global radial chromatin distribution to variations in the strength of TAD–TAD interactions compared to the lamin depleted nuclei, where a small (0.5 kT) increase of cross-type TAD–TAD interactions doubles the chromatin density in the central nucleus region. Conclusions A dynamical model of the entire fruit fly genome makes multiple genome-wide predictions of biological interest. The distribution of LADs along the chromatin chains affects their probabilities to be in contact with the NE and radial positioning of highly mobile TADs, playing a notable role in creating a non-random average global structure of the chromatin. We conjecture that an important role of attractive LAD–NE interactions is to stabilize global chromatin structure against inevitable cell-to-cell variations in TAD–TAD interactions. Supplementary Information The online version contains supplementary material available at 10.1186/s13072-023-00492-9.

It was suggested in [3] that the conformational state of chromatin across organisms depends on the chromosome to nucleus volume ratio, which emphasizes the need to make sure that computational models get this critical parameter right, in order to faithfully describe reality.

Bead-bead interactions and bead types
Bonded interactions between neighbouring beads along the chromatin chains are modeled by a harmonic potential, where r ij is the center-to-center distance between beads i and j = i + 1, the equilibrium bond length, l ij = r i + r j , is a sum of beads' radii, and κ is a bond strength parameter. Following Ref. [4], we choose κ = 10 kT to suppress chain crossing to a large degree. Non-bonded interactions between beads are described by a shortrange Lennard-Jones-cosine (LJ-cos) potential [5], with a minimum at r ij = l ij and σ ij = l ij /2 1 6 .
At r ij > l ij the LJ-cos function smoothly goes to 0 at 3 1 6 l ij , which is fulfilled by the corresponding choice of the parameters α and β (see Table 2 in Ref. [4]). Following our TAD-bead equivalence, we specify four corresponding bead types: Active, Null, PcG and HP1. The number of TADs in each epigenetic class and the corresponding number of beads in each bead type in our model are shown in Fig. S1. Non-bonded attractive interactions (Eq. 2) between beads of the same type are characterised by its own, bead-type specific interaction well depth parameter.
Attractive interactions (Eq. 2) between the TADs of different types -cross-type interactions -are not type-specific and are characterised by a single well depth parameter.
Repulsive interactions of the nucleolus with all beads, except HET beads, are described by a pure repulsive LJ potential function [5]: U r (r ij ) = 0 , r ij > l ij .
Here l ij = r i + r nuc is a contact distance between bead i and the nucleolus (j) of radius r nuc = 0.333 µm [1], ϵ r = 3 kT , and σ ij = l ij /2 1 6 . To localize HET beads near the nucleolus (see Ref. [1]), their interactions with the nucleolus are set to be strongly attractive, described by a standard LJ potential with ϵ i,nuc = 6 kT well depth and σ i,nuc = (r i + r nuc )/2 1 6 .
Interactions in specific pairs of remote loci Experimentally, 268 pairs of remote chromatin loci, with higher than expected contact probabilities, have been identified from the Hi-C data [6]. We map these specific loci pairs onto the pairs of TADs and introduce special interactions between the corresponding beads. The pairs of are ordered according to their degrees of contact enrichment: a ratio of the observed to expected contact probabilities, P obs ij /P exp ij . Since it is reasonable to assume that the bead-bead contact probability P ij ∼ exp(−∆U ij /kT ), where ∆U ij is a well depth of the effective interaction potential between beads i and j in the specific pair, one can express ln(P obs ij /P exp ij ) as Here, ∆U std is some unknown average "standard" attractive interaction between beads (TADs), assumed in the expected P exp ij . To estimate ∆U ij we fitted the experimental decay curve (see Fig. S2) of the logarithm of ordered (index n) degrees of enrichment for the 268 specific pairs of beads (TADs) with the formula: where index n of the ordered pairs runs from 0 to 267 and parameter C = 2.22.
Eqs. 4 and 5 allow us to determine the well depth ϵ n of the LJ potential used to describe attractive interactions in each specific pair (in kT units): Parameter E sp absorbs the constant C from Eq. 5 and the unknown interaction ∆U std from Eq. 4: E sp = C − ∆U std . We vary E sp between 3 and 6 kT and find that the optimum for the model TAD-TAD contact probability matrix (Hi-C map) ln (P_obs / P_exp) experimental data points approximating equation Figure S2 Determining the interaction strength between the special pairs of TADs. Shown is the logarithm of the ordered (index n) degree of enrichment (ln(P obs n /P exp n )) for the contact probabilities of 268 specific pairs of remote chromatin loci [6] and approximating fit (Eq. 5) used to determine the attractive interaction LJ well depth for the interactions between the beads in that specific pairs. corresponds to E sp = 3.5 kT . With this value of E sp , the attractive interaction well depth ϵ n varies from 2.3 to 3.0 kT for most (263 of 268) specific TAD pairs.

Nuclear envelope and its interactions with chromatin chains
We model the NE as a spherical boundary which confines the motion of the chromosomes and attracts LADs [7,8,9,10]. Three radii of the NE are used: 1.94 µm, 2.0 µm and 2.06 µm. These account for natural cell-to-cell variability; without it, chromatin density profiles near the NE are unrealistic compared to experiment. Following [4], we map positions of 412 LADs [8] (a median size of LADs is about 90 kb) onto the chains of 1169 TADs. If a TAD contains LAD, then the corresponding bead can attractively interact with the NE. After mapping, we have determined 350 TADs (L-TADs) that contain LADs. The fraction of L-TADs with two or more LADs is about 13%. The distribution of TADs among epigenetic classes and the numbers of L-TADs in each class in our model are shown in Fig. S1.
We describe all L-TAD-NE attractive interactions by the LJ-cos potential (Eq. 2) with a single well depth parameter ϵ L (neglecting the variations in the LAD-NE affinity due to different LAD sizes). The interaction distance r ij is defined as a distance r in from the bead i center to a spherical surface positioned at r n = 0.04 µm outside the NE boundary. The energy minimum is placed at the distance r i (bead i radius) from the NE, resulting in parameters l in = r i + r n and σ in = l in 2 1 6 .
Fraction of LADs at the NE, which depends on the affinity ϵ L , is calculated as the average (over the described above ensemble of 18 nuclei) fraction of L-TAD beads within (position of the bead centers) 0.09 µm layer (average bead radius) at the NE. The thickness of this layer corresponds to about 0.2 µm layer of beads which are in contact with the NE.
Interactions of the NE with beads not containing LADs are described by the purely repulsive potential U r (Eq. 3) with ϵ r = 3 kT , and with r ij , σ ij and l ij replaced by the L-TAD-NE interaction parameter r in , σ in and l in , defined above.

Lamin Mutant model
The Lamin mutant nucleus model has the L-TAD-NE affinity parameter ϵ L reduced to 0.1 kT level. The small but non-zero value is chosen to use the same LJ-cos function (Eq. 2) as for the other values of the L-TAD-NE affinity, and to reflect possible other tethering mechanisms of LADs to the NE [11] in the Lamin B (Dm0) knockdowns or mutants with much smaller levels of LAD-NE attraction than in the WT nuclei. We use Lamin B mutants as the experimental counterpart to our Lamin mutant nucleus model. All other interaction parameters are unchanged from the wild-type (WT) model.

Experimental chromatin density profiles
The results of the simulations are compared with the experimental data published previously in Ref. [12]. In that experiment, D. melanogaster line with the Lamin B (Dm0) mutant genotype Lam A25 pr 1 /CyO; st 1 (BSC, #25092) was used to investigate the role of the NE in nuclear architecture. Proventriculi (the end parts of the larval foreguts) were dissected from the Lam A25 mutant and Canton-S wild-type (WT) 3rd instar larvae. The radial distributions of DAPI-stained chromatin inside the nuclei were measured by quantification of the fluorescence intensity along the nuclei radii using confocal microscopy. The radial coordinates for each nucleus in pixels were translated into relative coordinates and normalized, see Ref. [12]. Mutation Lam A25 removes the CaaX box from Lamin B [13] disrupting the interaction of Lamin B with the NE. As a result, tethering of the chromatin to the nuclear periphery is impaired in the proventriculus nuclei of Lam A25 mutants as seen from the change in the chromatin radial distribution [12].

Matching simulation time to biological time
To relate the simulation timescale with the experiment [14], for the chromatin dynamics in D. melanogaster interphase nuclei, we compare [4] a diffusive motion of model beads in the simulations with the experimental interphase chromatin diffusion, observed in experiment [14], and use the match to estimate the scaling factor λ that converts the simulation time to real biological time. To characterize the diffusion, we calculate time dependencies of the distance R i (t) between a selected bead i and the nucleus center. Trajectories of for 9 randomly selected beads, which do not contain LADs and, thus, cannot be directly affected by the attraction to the NE, have been extracted. The dependence of the mean squared displacement (MSD) on the time interval ∆t, is calculated for 3 different values of the friction parameter γ (0.01/τ , 0.1/τ and 1/τ ). The averaging is performed over 9 selected beads along 18 trajectories for each value of γ. We fit the first 300 · 10 3 time steps of each of the three curves ⟨∆R 2 i (∆t)⟩ (see Fig. S3) with the equation for sub-diffusive motion of chromosomal loci [14], which describes MSD as a function of time interval ∆t, where D app is the apparent diffusion coefficient. This equation (with λ = 1 and ∆t in seconds) describes the experimentally observed diffusion of chromosomal loci over the time periods ranging from 1 to 10 3 s [14]. We obtain reasonable fits of the simulation derived MSDs with 4D app = 0.061 µm 2 and λ = 10 4 s −1 for γ = 0.01/τ , λ = 2.

Contact probability (Hi-C) map
The TAD-TAD contact probability map (Hi-C map) of the model is defined as the average of the 18 Hi-C maps (1169x1169) calculated over the 18 trajectories of the model nuclei of different sizes and mutual arrangements of the chromosomes. Each matrix element characterizes a contact probability p ij of the two beads defined as the fraction of the configurations where the beads i and j are in a contact. We define the bead-bead contact as the configuration with the distance between the centers of the two beads d ij < r i + r j + 0.2 µm.

Model development
Interaction parameters of the models Three stages of model development are used. In the initial stage, only HP1 type beads -classical heterochromatin TADs -were included into B-type set of beads. The rest of the beads were characterised by a weak 0.1 kT TAD-TAD attractive interaction. We analyzed 144 combinations of the interaction parameters by generating 40 × 10 6 time step MD trajectories for each parameter set, using only one nucleus topology. We scanned the following range of parameters with in 1 kT increments. The attractive LJ-cos interactions: 3 -6 kT for the E sp parameter in the specific TAD pair interactions; 1 -6 kT for the HP1-HP1 (B-type beads) interactions; 2 -6 kT for the LAD-NE (L-TAD-NE) interaction parameter ϵ L . For each parameter set the Hi-C map was calculated and compared with the experimental Hi-C map [6,1] resulting in the range of Pearson's correlation coefficients from 0.895 to 0.950. This analysis allowed us to narrow down the ranges of the parameters for further refinement.
In the next stage of our model development, in addition to the HP1 beads, we considered transcriptionally inactive [6] Null and PcG beads as members of the inactive chromatin B-compartments. The corresponding attractive interactions (Eq. 2) between beads of the same type in these B-compartments are characterized by the LJ-cos well depth parameters ϵ B . To allow possible compartmentalization within each of these three bead types, we describe the cross-type interactions between these beads (Eq. 2) by the same small LJ-cos parameter as the one, ϵ AB , used for the A-B cross-type interactions. Active-Active interactions (Eq. 2) between beads in A-type compartments are described by the LJ-cos parameter ϵ A . Interactions between HET and CEN beads are considered as A-B cross-type interactions.
Results from the previous stage allowed us to focus on a more narrow range of model parameters: only 47 sets of the interaction parameters. The selected ranges for the LJ-cos interaction well depths are: 1 -2 kT for ϵ B ; 3 -4 kT for E sp ; and 2 -6 kT for the ϵ L (L-TAD-NE interactions). Parameters ϵ A and ϵ AB were fixed at 0.1 kT (to satisfy selection rule #3 above for ϵ AB ). For each set, a 40×10 6 time-step trajectory was generated for a single CIS-X6S nucleus topology (see Fig. 1 in "Methods", main text), the Hi-C map was calculated and compared with the experimental Hi-C map [6,1]. The comparisons produced a range of Pearson's correlation coefficients from 0.932 to 0.953. The radial distributions of L-TADs were calculated and the fractions of L-TADs at the NE were compared with the experimental data [15] (selection rule #2 in main text, "Model Development").
To further refine the interaction parameters, in the last stage of model development, we performed simulations of the full ensembles of nuclei topologies (see Fig. 1 in the main text), properly weighted as described in "Methods" (main text), including the use of 3 different (2.00 ± 0.06 µm) nucleus sizes. As described in "Methods" (main text), the resulting ensemble of simulated nuclei consists of 18 different initial states and the corresponding trajectories; the latter are used to obtain the topology and size averaged chromatin distributions and Hi-C maps for each test combination of the interaction parameters. At this stage, 32 parameter sets have been used to determine the "optimum" WT parameter set and to investigate the effects of parameter deviations onto the chromatin properties. 40 × 10 6 time-steps trajectories were generated for each of characterised by a single parameter set. Compared to the previous optimization stage, and based on its results, a more narrow range (3 -4.5 kT ) of the LAD-NE interactions parameters was used. Several values of the cross-type ϵ AB parameter and ϵ A parameter (in the range 0.1 -1.0 kT ) were used as well. In addition, 8 sets of the interaction parameters for the Lamin mutant model were considered to investigate its chromatin properties.
After calculating the ensemble averaged Hi-C maps for each parameter set, we selected 9 WT-like sets with the highest (greater than 0.955) Pearson's correlation coefficients with the experimental Hi-C map (selection rule #1 in the main text). Given the range (1 -2 kT ) used for the interaction parameter ϵ B , the selection rule #3 is satisfied when the A-B cross-type interaction parameter ϵ AB does not exceed 0.5 kT . This reduces the number of the selected sets from 9 to 5. All these new 5 sets have the same ϵ AB = 0.5 kT and all 3 types of ϵ B parameters (for Null, PcG and HP1 beads) equal to 1.5 kT , and are characterized by roughly the same Pearson's correlation coefficients, ∼ 0.956. However, the values of the LAD-NE interaction parameters ϵ L are still different (3.0, 3.5 and 4.0 kT ), allowing us to apply the selection rule #2 for further refinement. Slightly different are the values of the ϵ A parameter (0.1 and 0.5 kT ) and E sp parameter (3.0 and 3.5 kT ), suggesting that these differences don't have a noticeable effect on the local chromatin structure.
We would like to note here that the Pearson's correlation coefficients for the ensemble averaged model Hi-C maps for these "best" 5 parameter sets are higher than the coefficients for each individual nucleus Hi-C map in the ensemble of nuclei. Most of these single nucleus Hi-C coefficients do not exceed the largest value obtained in the previous parameter optimization stage.
Selecting the A-A interaction parameter ϵ A to be 0.1 kT , as found in the 4 of 5 sets, and having it lower than ϵ AB = 0.5 kT for the cross-type interactions between different classes of B-type beads (e.g., HP1-PcG interactions), and excluding the only set with E sp = 3.0 kT (which has the smallest Pearson's correlation coefficient), we have determined a single set of the parameters (other than ϵ L ) we can use to select the LAD-NE interaction strength ϵ L in our WT nucleus model.
To select a single value of the ϵ L parameter (ϵ = ϵ L in Eq. 2) that reproduces experimentally observed fraction of LADs near the NE (25% [15], selection rule #2), we analyze the radial distributions of L-TADs in the nuclei models with the selected above set of the interaction parameters. The normalized cumulative distributions shown in Fig. S4 suggest that the optimum value of the L-TAD-NE interaction well depth ϵ L is 4.0 kT . This value ensures that, on average, 25% of L-TAD bead centers are within 0.09 µm (average bead radius) from the average position of the NE (2.0 µm from the nucleus center). Experimental data [15] show that about 25% (horizontal dashed line at 0.75 mark) of LADs are on average located at the NE.
The final optimum WT set of the interaction parameters is shown in the second column of the Table 1 in the main text. We use this set for 10x longer simulations reported in the next section. The resulting chromatin density profiles will be compared with the available experimental data [12]. We will also investigate how sensitive is the chromatin structure to the deviations in the parameters.

Compartmentalization
To see whether our model produces the expected different levels of compartmentalization for A and B types of beads, we used the model derived Hi-C and estimated contact intensities for two types of beads -Active beads (A type) and Null beads (B type). For estimation of the contact intensities we calculated average sums of the contact probabilities for these two types of beads with all other surrounding beads.
The resulting values of the sums -3.86 for Active beads and 5.83 for Null beadssuggest that compartmentalization, characterized by the intensities of contacts, for A type beads is roughly twice less intense than for B type beads.
Heterogeneity of the biological data used in the model Drosophila nuclei at interphase are modeled using a heterogeneous set of experimental data, which were available to us. The Hi-C data are from late Drosophila embryos [6]. The Lamin B mutant data are from larval low-level polytene nuclei that change chromatin distribution similarly to regular diploid nuclei [12,16]. The LAD data are from Drosophila Kc cells of embryonic origin [15,8]. Since these data are not from highly-specialized tissues, we believe they represent a typical Drosophila nucleus. Moreover, LADs appear to be robust in different D. melanogaster cells including polytene chromosomes [15].

Results
Temporal evolution of the WT model Hi-C maps NE as an "attractive enclosure" The chromatin distribution in the Lamin mutant nuclei resembles a more condensed, globule-like form compared to WT nuclei. A question arises: does the NE in WT nuclei acts only as an attractive enclosure needed to "keep interphase chromosomes slightly stretched" [16] on biologically relevant time-scales? Will the chromatin condense in the complete absence of the NE, not only in the absence of LAD-NE affinity? To answer these questions we have carried out a series of "what if" simulations of the model nuclei with the NE completely removed. A typical result is shown in Fig. S6. The chromosomes partially decondense and extend substantially, and the chromosome chains partially detach from each other resembling an "extended coil" polymer state, on the time scale of a minute. Thus, the NE, even without attractive interactions with LADs, is still necessary to suppress the fluctuations of the chromatin preventing it from unfolding, very quickly. Figure S6 The effect of complete removal of the NE: chromatin decompacts. Within 1 minute, chromosome 4 (cyan) dissociates and drifts away from the rest three de-condensed chromosomes.

Figure S8
Difference Hi-C map between Lamin mutant model Hi-C map (Fig. 5, main text) and the corresponding WT model Hi-C map (Fig. 3, main text, bottom panel). Note that the intensity scale is 10 times smaller here than that of the original Hi-C maps, pointing to relatively small differences in TAD-TAD contact probabilities. Most TAD-TAD contacts, including inter-chromosome contacts, are slightly enhanced (red areas) in the Lamin mutant compared to the WT nuclei. Small areas where the contact frequency decreases (blue spots) in the Lamin mutant are limited to close to diagonal intra-arm contacts.  (cytological region 60D in Ref. [16]) in the WT and Lamin mutant nucleus models.